This document contains all code necessary to replicate the analyses for the above article. Any questions about this project can be directed to Kaysee Arrowsmith at kcarrows@utexas.edu.
library(knitr)
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
library(tidyverse)
theme_set(theme_classic())
library(vegan)
library(ecodist)
library(mclogit)
library(glmmTMB)
library(DHARMa)
library(RColorBrewer)
library(geosphere)
library(Matrix)
library(lme4)
library(lmerTest)
library(ggeffects)
MRM.base <- function (formula = formula(data), data, nperm = 1000, method = "linear", mrank = FALSE, permute.control = "permute.control") {
m <- match.call(expand.dots = FALSE)
m2 <- match(c("formula", "data"), names(m), nomatch = 0)
m <- m[c(1, m2)]
m[[1]] <- as.name("model.frame")
m <- eval(m, parent.frame())
m <- as.matrix(m)
n <- (1 + sqrt(1 + 8 * nrow(m)))/2
if (abs(n - round(n)) > 1e-07)
stop("Matrix not square.\n")
n <- round(n)
if (ncol(m) < 2)
stop("Not enough data. \n")
if (method == "linear") {
X <- m[, 2:ncol(m), drop = FALSE]
Y <- m[, 1, drop = FALSE]
colnames(Y) <- "Y"
newdata <- data.frame(Y = Y, X)
fit1 <- lm(Y ~ ., data = newdata)
b <- coefficients(fit1)
sig <- summary(fit1)$sigma
sig.df <- summary(fit1)$df[2]
b.pval <- NA
dev.pval <- NA
if (nperm > 0) {
b.all <- matrix(NA, nrow = nperm, ncol = length(b))
b.all[1, ] <- b
sig.all <- rep(NA, nperm)
sig.all[1] <- sig
for (i in 2:nperm) {
newSample <- shuffle(n, control = permute.control)
newY <- full(Y)
newY <- newY[newSample, newSample]
newY <- lower(newY)
newdata <- data.frame(Y = newY, X = X)
newfit <- lm(Y ~ ., data = newdata)
b.all[i, ] <- coefficients(newfit)
sig.all[i] <- summary(newfit)$sigma
}
b.pval <- apply(b.all, 2, function(x) length(x[abs(x) >=
abs(x[1])])/nperm)
R2.pval <- length(sig.all[sig.all >= sig.all[1]])/nperm
}
results <- list(coef = cbind(b, pval = b.pval),
R2 = c(R2 = sig, R2.df = sig.df,
R2.pval = R2.pval))
}
else {
if (method == "logistic") {
X <- m[, 2:ncol(m), drop = FALSE]
Y <- m[, 1, drop = FALSE]
colnames(Y) <- "Y"
newdata <- data.frame(Y = Y, X)
fit1 <- glm(Y ~ ., data = newdata, family = binomial(link = "logit"))
b <- coefficients(fit1)
dev <- summary(fit1)$deviance
dev.df <- summary(fit1)$df.residual
b.pval <- NA
dev.pval <- NA
if (nperm > 0) {
b.all <- matrix(NA, nrow = nperm, ncol = length(b))
b.all[1, ] <- b
dev.all <- rep(NA, nperm)
dev.all[1] <- dev
for (i in 2:nperm) {
newSample <- shuffle(n, control = permute.control)
newY <- full(Y)
newY <- newY[newSample, newSample]
newY <- lower(newY)
newdata <- data.frame(Y = newY, X = X)
newfit <- glm(Y ~ ., data = newdata, family = binomial(link = "logit"))
b.all[i, ] <- coefficients(newfit)
dev.all[i] <- summary(newfit)$deviance
}
b.pval <- apply(b.all, 2, function(x) length(x[abs(x) >=
abs(x[1])])/nperm)
dev.pval <- length(dev.all[dev.all >= dev.all[1]])/nperm
}
results <- list(coef = cbind(b, pval = b.pval),
dev = c(resid.dev = dev, resid.df = dev.df,
dev.pval = dev.pval))
}
else {
stop("method must be 'linear' or 'logistic'\n")
}
}
results
}
dat <- read.csv("data/tempint-networks.csv", stringsAsFactors = F) %>%
mutate(date = as.Date(date, "%F"),
year = format(date, "%Y"))
plants <- read.csv("data/tempint-plants.csv", stringsAsFactors = F) %>%
mutate(date = as.Date(date, "%F"),
year = format(date, "%Y"))
samptemps <- read.csv("data/tempint-samptemps.csv", stringsAsFactors = F) %>%
mutate(date = as.Date(date, "%F"),
year = format(date, "%Y")) %>%
filter(!is.na(round),
site %in% dat$site)
siteinfo <- read.csv("data/tempint-siteinfo.csv", stringsAsFactors = F)
Calculate every possible interaction that could have occurred based on the plants and pollinators observed at that date/site/round.
# Family
poss.family <- dat %>%
filter(insect_family != "") %>%
group_by(date, site, round, insect_family) %>%
tally() %>%
full_join(plants, by = c("date", "site"), multiple = "all") %>%
dplyr::select(date, site, round, insect_family, plant) %>%
mutate(concat = as.factor(paste(date, site, round, sep = ".")),
date.site = as.factor(paste(date, site, sep = ".")))
# Genus
poss.genus <- dat %>%
filter(insect_genus != "") %>%
group_by(date, site, round, insect_genus) %>%
tally() %>%
full_join(plants, by = c("date", "site"), multiple = "all") %>%
dplyr::select(date, site, round, insect_genus, plant) %>%
mutate(concat = as.factor(paste(date, site, round, sep = ".")),
date.site = as.factor(paste(date, site, sep = ".")))
# Functional Group
poss.func <- dat %>%
filter(gross_ID != "Unknown") %>%
mutate_if(grepl("Bombus*", .), ~replace(., grepl("Bombus*", .), "Bombus")) %>%
group_by(date, site, round, gross_ID) %>%
tally() %>%
full_join(plants, by = c("date", "site"), multiple = "all") %>%
dplyr::select(date, site, round, gross_ID, plant) %>%
mutate(concat = as.factor(paste(date, site, round, sep = ".")),
date.site = as.factor(paste(date, site, sep = ".")))
Which of the possible interactions (above) actually occurred.
# Family
occ.family <- dat %>%
filter(insect_family != "") %>%
group_by(date, site, round, insect_family, plant) %>%
tally(name = "int") %>%
right_join(poss.family, by = c("date", "site", "round", "insect_family", "plant")) %>%
mutate(int = ifelse(is.na(int), 0, int))
# Genus
occ.genus <- dat %>%
filter(insect_genus != "") %>%
group_by(date, site, round, insect_genus, plant) %>%
tally(name = "int") %>%
right_join(poss.genus, by = c("date", "site", "round", "insect_genus", "plant")) %>%
mutate(int = ifelse(is.na(int), 0, int))
# Functional Group
occ.func <- dat %>%
filter(gross_ID != "Unknown") %>%
mutate_if(grepl("Bombus*", .), ~replace(., grepl("Bombus*", .), "Bombus")) %>%
group_by(date, site, round, gross_ID, plant) %>%
tally(name = "int") %>%
right_join(poss.func, by = c("date", "site", "round", "gross_ID", "plant")) %>%
mutate(int = ifelse(is.na(int), 0, int))
Attach the temperature for each sampling round. We also calculate some relative abundances here that will be useful for plotting our results later on.
# Family
final.family <- occ.family %>%
left_join(samptemps, by = c("date", "site", "round")) %>%
mutate(int.ID = paste(insect_family, plant, sep = ".")) %>%
dplyr::select(date, site, round, date.site, concat, insect_family, plant, int.ID, int, avg.temp) #, latitude, longitude)
# Genus
final.genus <- occ.genus %>%
left_join(samptemps, by = c("date", "site", "round")) %>%
mutate(int.ID = paste(insect_genus, plant, sep = "."),
year = format(date, "%Y")) %>%
# midtime = as.numeric(strptime(starttime, format = "%H:%M")) + as.numeric(strptime(endtime, format = "%H:%M"))/2,
# doy = format(date, "%j")) %>%
dplyr::select(date, year, site, round, date.site, concat, insect_genus, plant, int.ID, int, avg.temp)#, latitude, longitude, midtime, doy)
final.func <- occ.func %>%
left_join(samptemps, by = c("date", "site", "round")) %>%
mutate(int.ID = paste(gross_ID, plant, sep = "."),
year = format(date, "%Y")) %>%
# midtime = as.numeric(strptime(starttime, format = "%H:%M")) + as.numeric(strptime(endtime, format = "%H:%M"))/2,
# doy = format(date, "%j")) %>%
dplyr::select(date, year, site, round, date.site, concat, gross_ID, plant, int.ID, int, avg.temp)#, latitude, longitude, midtime, doy)
This checks whether temperature varied with round. A priori, we would assume that temperature varies with date, but that should not matter in our analysis since our choice sets are always within a single date. However, we include date here too in order to control for that effect.
round.mod <- lmer(avg.temp ~ date + as.factor(round) + (1|site), data = samptemps)
summary(round.mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg.temp ~ date + as.factor(round) + (1 | site)
## Data: samptemps
##
## REML criterion at convergence: 1156.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.66524 -0.58990 -0.04593 0.65944 2.33962
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.275 1.508
## Residual 16.298 4.037
## Number of obs: 203, groups: site, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 101.689518 31.697145 197.176783 3.208 0.00156 **
## date -0.004049 0.001663 197.414224 -2.436 0.01576 *
## as.factor(round)2 2.597828 0.680227 194.334156 3.819 0.00018 ***
## as.factor(round)3 3.253508 0.704645 194.432966 4.617 7.06e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) date as.()2
## date -1.000
## as.fctr(r)2 -0.015 0.005
## as.fctr(r)3 -0.022 0.012 0.487
anova(round.mod)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## date 96.68 96.684 1 197.41 5.9321 0.01576 *
## as.factor(round) 400.27 200.134 2 194.37 12.2795 9.508e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
samptemps.BC21 <- samptemps %>% filter(site == "Brush Creek Hill", year == "2021")
test.bc21 <- aov(avg.temp ~ as.factor(round), data = samptemps.BC21)
summary(test.bc21)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 6.32 3.158 0.18 0.837
## Residuals 14 245.65 17.547
samptemps.BC22 <- samptemps %>% filter(site == "Brush Creek Hill", year == "2022")
test.bc22 <- aov(avg.temp ~ as.factor(round), data = samptemps.BC22)
summary(test.bc22)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 61.9 30.93 1.184 0.333
## Residuals 15 391.7 26.11
samptemps.KP21 <- samptemps %>% filter(site == "Kebler Pass", year == "2021")
test.kp21 <- aov(avg.temp ~ as.factor(round), data = samptemps.KP21)
summary(test.kp21)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 78.25 39.12 2.197 0.142
## Residuals 17 302.75 17.81
samptemps.KP22 <- samptemps %>% filter(site == "Kebler Pass", year == "2022")
test.kp22 <- aov(avg.temp ~ as.factor(round), data = samptemps.KP22)
summary(test.kp22)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 99.6 49.78 1.409 0.271
## Residuals 17 600.4 35.32
samptemps.OPH22 <- samptemps %>% filter(site == "Ohio Pass High", year == "2022")
test.oph22 <- aov(avg.temp ~ as.factor(round), data = samptemps.OPH22)
summary(test.oph22)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 66.87 33.44 4.459 0.0277 *
## Residuals 17 127.47 7.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
samptemps.OPL22 <- samptemps %>% filter(site == "Ohio Pass Low", year == "2022")
test.opl22 <- aov(avg.temp ~ as.factor(round), data = samptemps.OPL22)
summary(test.opl22)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 162.2 81.11 6.06 0.00972 **
## Residuals 18 240.9 13.39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
samptemps.RG21 <- samptemps %>% filter(site == "Rustlers Gulch", year == "2021")
test.rg21 <- aov(avg.temp ~ as.factor(round), data = samptemps.RG21)
summary(test.rg21)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 23.24 11.621 1.563 0.237
## Residuals 18 133.82 7.434
samptemps.RG22 <- samptemps %>% filter(site == "Rustlers Gulch", year == "2022")
test.rg22 <- aov(avg.temp ~ as.factor(round), data = samptemps.RG22)
summary(test.rg22)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 26.69 13.343 2.097 0.149
## Residuals 20 127.25 6.363
samptemps.SF21 <- samptemps %>% filter(site == "Stupid Falls", year == "2021")
test.sf21 <- aov(avg.temp ~ as.factor(round), data = samptemps.SF21)
summary(test.sf21)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 14.04 7.02 0.435 0.654
## Residuals 17 274.08 16.12
samptemps.SF22 <- samptemps %>% filter(site == "Stupid Falls", year == "2022")
test.sf22 <- aov(avg.temp ~ as.factor(round), data = samptemps.SF22)
summary(test.sf22)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round) 2 31.71 15.85 1.141 0.339
## Residuals 20 277.90 13.89
Round does significantly impact temperature, so it will make sense to include round as a main effect in our subsequent analyses.
We use a mixed conditional logit model (mclogit) to
determine whether pollinators’ foraging preferences are affected by
temperature. In this model, we can imagine each interaction as a
discrete choice made by a pollinator among the choice set of all flowers
available at that site/date combination (defined by
site.date). We include round as a random intercept and the
identities of the insect and plant involved in the interaction as random
slopes.
First, we run this analysis with insects grouped by functional group. This comprises the largest amount of data.
mc.func <- mclogit(cbind(int, as.integer(date.site)) ~
poly(avg.temp, 2, raw = TRUE) +
as.factor(round),
random = list(~1 + avg.temp|gross_ID,
~1 + avg.temp|plant),
data = final.func)
##
## Iteration 1 - deviance = 7609.395 - criterion = 0.9746249
## Iteration 2 - deviance = 7376.174 - criterion = 0.02211211
## Iteration 3 - deviance = 7343.175 - criterion = 0.00640909
## Iteration 4 - deviance = 7331.056 - criterion = 0.002133236
## Iteration 5 - deviance = 7327.29 - criterion = 0.0002740379
## Iteration 6 - deviance = 7326.449 - criterion = 6.167405e-06
## Iteration 7 - deviance = 7326.374 - criterion = 2.948911e-09
## converged
summary(mc.func)
##
## Call:
## mclogit(formula = cbind(int, as.integer(date.site)) ~ poly(avg.temp,
## 2, raw = TRUE) + as.factor(round), data = final.func, random = list(~1 +
## avg.temp | gross_ID, ~1 + avg.temp | plant))
##
## Coefficents:
## Estimate Std. Error z value Pr(>|z|)
## poly(avg.temp, 2, raw = TRUE)1 0.224690 0.073511 3.057 0.00224 **
## poly(avg.temp, 2, raw = TRUE)2 -0.005419 0.001262 -4.294 1.76e-05 ***
## as.factor(round)2 0.106515 0.060862 1.750 0.08010 .
## as.factor(round)3 -0.166317 0.070442 -2.361 0.01822 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Co-)Variances:
## Grouping level: gross_ID
## Estimate Std.Err.
## (Const.) 3.926835 2.935e+00
## avg.temp -0.005289 0.004275 3.952e-03 5.306e-06
##
## Grouping level: plant
## Estimate Std.Err.
## (Const.) 2.186126 1.03431
## avg.temp 0.033240 0.004818 0.01576 0.00024
##
## Approximate residual deviance: 7326
## Number of Fisher scoring iterations: 7
## Number of observations
## Groups by gross_ID: 15
## Groups by plant: 64
## Individual observations: 2585
Next, we run this analysis with insects identified to the family level. This is a coarser-scale analysis than genus level, but there are a larger number of data points at the family level (because it includes insects for which we can identify to family but not to genus).
mc.family <- mclogit(cbind(int, as.integer(date.site)) ~
poly(avg.temp, 2, raw = TRUE) +
as.factor(round),
random = list(~1 + avg.temp|insect_family,
~1 + avg.temp|plant),
data = final.family)
##
## Iteration 1 - deviance = 6469.196 - criterion = 0.9729606
## Iteration 2 - deviance = 6242.679 - criterion = 0.110585
## Iteration 3 - deviance = 6210.56 - criterion = 0.009696167
## Iteration 4 - deviance = 6198.641 - criterion = 0.002914789
## Iteration 5 - deviance = 6194.957 - criterion = 0.0003635539
## Iteration 6 - deviance = 6194.189 - criterion = 7.550973e-06
## Iteration 7 - deviance = 6194.122 - criterion = 2.866196e-09
## converged
summary(mc.family)
##
## Call:
## mclogit(formula = cbind(int, as.integer(date.site)) ~ poly(avg.temp,
## 2, raw = TRUE) + as.factor(round), data = final.family, random = list(~1 +
## avg.temp | insect_family, ~1 + avg.temp | plant))
##
## Coefficents:
## Estimate Std. Error z value Pr(>|z|)
## poly(avg.temp, 2, raw = TRUE)1 0.170662 0.079078 2.158 0.030916 *
## poly(avg.temp, 2, raw = TRUE)2 -0.004732 0.001390 -3.404 0.000664 ***
## as.factor(round)2 0.104609 0.069255 1.510 0.130917
## as.factor(round)3 -0.106748 0.080308 -1.329 0.183770
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Co-)Variances:
## Grouping level: insect_family
## Estimate Std.Err.
## (Const.) 1.88751 1.6661393
## avg.temp 0.03328 0.00509 0.0294493 0.0005205
##
## Grouping level: plant
## Estimate Std.Err.
## (Const.) 1.971802 1.119256
## avg.temp 0.038663 0.004551 0.021988 0.000432
##
## Approximate residual deviance: 6194
## Number of Fisher scoring iterations: 7
## Number of observations
## Groups by insect_family: 37
## Groups by plant: 64
## Individual observations: 2124
We find that temperature has highly significant linear and quadratic effects.
Next, we run this analysis again but at the genus level. This is a finer-scale analysis that eliminates some uncertainty associated with the potential for variation among genera within the same family, but it loses some statistical power since we have roughly 500 fewer interactions that can be identified to the genus level.
mc.genus <- mclogit(cbind(int, as.integer(date.site)) ~
poly(avg.temp, 2, raw = TRUE) +
round,
random = list(~1 + avg.temp|insect_genus,
~1 + avg.temp|plant),
data = final.genus)
##
## Iteration 1 - deviance = 5189.292 - criterion = 0.9178357
## Iteration 2 - deviance = 5086.866 - criterion = 0.02430653
## Iteration 3 - deviance = 5058.967 - criterion = 0.003680271
## Iteration 4 - deviance = 5029.556 - criterion = 0.001792912
## Iteration 5 - deviance = 5000.943 - criterion = 0.001973756
## Iteration 6 - deviance = 4974.1 - criterion = 0.001997291
## Iteration 7 - deviance = 4950.184 - criterion = 0.001860424
## Iteration 8 - deviance = 4929.955 - criterion = 0.001618803
## Iteration 9 - deviance = 4913.621 - criterion = 0.001343169
## Iteration 10 - deviance = 4900.911 - criterion = 0.001084788
## Iteration 11 - deviance = 4891.257 - criterion = 0.0008686005
## Iteration 12 - deviance = 4884.002 - criterion = 0.0007008919
## Iteration 13 - deviance = 4878.535 - criterion = 0.0005777871
## Iteration 14 - deviance = 4874.354 - criterion = 0.0004907416
## Iteration 15 - deviance = 4871.082 - criterion = 0.0004296845
## Iteration 16 - deviance = 4868.451 - criterion = 0.0003847052
## Iteration 17 - deviance = 4866.282 - criterion = 0.0003468096
## Iteration 18 - deviance = 4864.461 - criterion = 0.0003081801
## Iteration 19 - deviance = 4862.92 - criterion = 0.0002624472
## Iteration 20 - deviance = 4861.631 - criterion = 0.0002056703
## Iteration 21 - deviance = 4860.59 - criterion = 0.0001392979
## Iteration 22 - deviance = 4859.925 - criterion = 5.395226e-05
## Iteration 23 - deviance = 4859.545 - criterion = 2.051158e-05
## Iteration 24 - deviance = 4859.299 - criterion = 2.425215e-05
## Iteration 25 - deviance = 4859.164 - criterion = 3.435256e-05
summary(mc.genus)
##
## Call:
## mclogit(formula = cbind(int, as.integer(date.site)) ~ poly(avg.temp,
## 2, raw = TRUE) + round, data = final.genus, random = list(~1 +
## avg.temp | insect_genus, ~1 + avg.temp | plant))
##
## Coefficents:
## Estimate Std. Error z value Pr(>|z|)
## poly(avg.temp, 2, raw = TRUE)1 0.162469 0.092899 1.749 0.0803 .
## poly(avg.temp, 2, raw = TRUE)2 -0.004019 0.001651 -2.435 0.0149 *
## round -0.029829 0.043986 -0.678 0.4977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Co-)Variances:
## Grouping level: insect_genus
## Estimate Std.Err.
## (Const.) 3.026e-02 3.860e-06
## avg.temp -1.238e-03 9.927e-05 1.582e-07 6.481e-09
##
## Grouping level: plant
## Estimate Std.Err.
## (Const.) 4.41538 1.0734239
## avg.temp -0.08288 0.00359 0.0201589 0.0003786
##
## Approximate residual deviance: 4859
## Number of Fisher scoring iterations: 25
## Number of observations
## Groups by insect_genus: 52
## Groups by plant: 64
## Individual observations: 1627
## Note: Algorithm did not converge.
We still see a highly significant quadratic effect of temperature, but only marginal significance in the linear effect of temperature. This loss of significance seems likely to be related to the loss of statistical power in this genus-level analysis.
Since we found similar results at both the family and genus levels, we perform all subsequent analyses at the genus level to take advantage of the finer taxonomic resolution.
We will check whether temperature affects the number of pollinators actively foraging during our sampling periods.
func.active <- final.func %>%
mutate(round = as.factor(round)) %>%
group_by(date, site, round, avg.temp, gross_ID) %>%
summarise(num.active = sum(int))
total.func <- func.active %>%
group_by(date, site, round) %>%
summarise(total.active = sum(num.active))
func.subset <- func.active %>%
group_by(date, site, gross_ID) %>%
tally() %>%
filter(n == 3) %>%
inner_join(func.active, by = c("date", "site", "gross_ID")) %>%
dplyr::select(-n) %>%
left_join(total.func, by = c("date", "site", "round"))
activity.func <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|gross_ID),
data = func.subset,
family = binomial)
summary(activity.func)
## Family: binomial ( logit )
## Formula:
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,
## 2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |
## gross_ID)
## Data: func.subset
##
## AIC BIC logLik deviance df.resid
## 1888.6 1927.6 -934.3 1868.6 353
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## date (Intercept) 1.090e-01 3.302e-01
## site (Intercept) 3.265e-09 5.714e-05
## gross_ID (Intercept) 4.904e+00 2.215e+00
## avg.temp 4.341e-03 6.589e-02 -0.95
## Number of obs: 363, groups: date, 46; site, 6; gross_ID, 7
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9604050 1.6569274 -0.580 0.562
## poly(avg.temp, 2, raw = TRUE)1 -0.0229519 0.1025394 -0.224 0.823
## poly(avg.temp, 2, raw = TRUE)2 0.0005544 0.0017317 0.320 0.749
## round2 -0.0753297 0.0815248 -0.924 0.355
## round3 -0.0773739 0.0887164 -0.872 0.383
actfunc.resid <- simulateResiduals(activity.func, n = 250)
plot(actfunc.resid)
We see no significance of any fixed effects. Diagnostic plots all look good.
family.active <- final.family %>%
mutate(round = as.factor(round)) %>%
group_by(date, site, round, avg.temp, insect_family) %>%
summarise(num.active = sum(int))
total.family <- family.active %>%
group_by(date, site, round) %>%
summarise(total.active = sum(num.active))
family.subset <- family.active %>%
group_by(date, site, insect_family) %>%
tally() %>%
filter(n == 3) %>%
inner_join(family.active, by = c("date", "site", "insect_family")) %>%
dplyr::select(-n) %>%
left_join(total.family, by = c("date", "site", "round"))
activity.family <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|insect_family),
data = family.subset,
family = binomial)
summary(activity.family)
## Family: binomial ( logit )
## Formula:
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,
## 2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |
## insect_family)
## Data: family.subset
##
## AIC BIC logLik deviance df.resid
## 1176.5 1211.6 -578.3 1156.5 236
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## date (Intercept) 2.168e-01 4.657e-01
## site (Intercept) 3.246e-09 5.698e-05
## insect_family (Intercept) 7.700e+00 2.775e+00
## avg.temp 1.149e-02 1.072e-01 -0.98
## Number of obs: 246, groups: date, 40; site, 6; insect_family, 14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3495175 2.1520437 -0.162 0.871
## poly(avg.temp, 2, raw = TRUE)1 -0.0362095 0.1420572 -0.255 0.799
## poly(avg.temp, 2, raw = TRUE)2 0.0002365 0.0024072 0.098 0.922
## round2 -0.1470876 0.1092397 -1.346 0.178
## round3 0.0963853 0.1209730 0.797 0.426
actfam.resid <- simulateResiduals(activity.family, n = 250)
plot(actfam.resid)
No significance of any fixed effects at the family-level. KS test is
showing significance consistent with overdispersion. We do know that our
data are overdispersed by the model will not converge with
family = betabinomial. Overdispersion should artificially
lower statistical significance (Type I error), so since we don’t see any
statistical significance, we think it’s fine.
genus.active <- final.genus %>%
mutate(round = as.factor(round)) %>%
group_by(date, site, round, avg.temp, insect_genus) %>%
summarise(num.active = sum(int))
total.genus <- genus.active %>%
group_by(date, site, round) %>%
summarise(total.active = sum(num.active))
genus.subset <- genus.active %>%
group_by(date, site, insect_genus) %>%
tally() %>%
filter(n == 3) %>%
inner_join(genus.active, by = c("date", "site", "insect_genus")) %>%
dplyr::select(-n) %>%
left_join(total.genus, by = c("date", "site", "round"))
activity.genus <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|insect_genus),
data = genus.subset,
family = binomial)
summary(activity.genus)
## Family: binomial ( logit )
## Formula:
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,
## 2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |
## insect_genus)
## Data: genus.subset
##
## AIC BIC logLik deviance df.resid
## 721.1 751.6 -350.5 701.1 146
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## date (Intercept) 0.67368 0.8208
## site (Intercept) 0.14004 0.3742
## insect_genus (Intercept) 19.01941 4.3611
## avg.temp 0.02255 0.1502 -0.97
## Number of obs: 156, groups: date, 33; site, 6; insect_genus, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.229362 3.385593 0.954 0.340
## poly(avg.temp, 2, raw = TRUE)1 -0.304036 0.223456 -1.361 0.174
## poly(avg.temp, 2, raw = TRUE)2 0.005317 0.003841 1.384 0.166
## round2 -0.049612 0.164785 -0.301 0.763
## round3 0.230728 0.190392 1.212 0.226
actgen.resid <- simulateResiduals(activity.genus, n = 250)
plot(actgen.resid)
No significance of any fixed effects at the genus level. As with the family-level diagnostics, we see significance in the KS test consistent with overdispersion. The genus-level analysis will converge with a betabinomial distribution but lack of significance remains the same.
Next, we will use matrix regression to examine how dissimilarity in community composition relates to dissimilarity in temperature. For this analysis, we will include site/year and round as controlling factors that limit how permutations are conducted. We do not include round as a fixed effect because the MRM does not work with categorical variables.
plants.sum <- plants %>%
filter(plant != "") %>%
mutate(site.date = paste(site, date, sep = ".")) %>%
group_by(site.date, plant) %>%
summarise(total_flowers = sum(total_flowers))
func.polmatrix <- final.func %>%
group_by(concat, gross_ID) %>%
summarise(num.active = sum(int)) %>%
pivot_wider(names_from = gross_ID,
values_from = num.active,
values_fill = 0) %>%
arrange(concat) %>%
column_to_rownames(var = "concat")
func.polbeta <- as.matrix(vegdist(func.polmatrix, method = "horn"))
func.explanatory <- data.frame(concat = rownames(func.polmatrix)) %>%
separate(concat, into = c("date", "site", "round"), sep = "[.]", remove = F) %>%
mutate(date = as.Date(date),
round = as.numeric(round),
site.date = paste(site, date, sep = ".")) %>%
left_join(samptemps, by = c("date", "site", "round")) %>%
mutate(doy = strftime(date, format = "%j"),
year = strftime(date, format = "%Y")) %>%
arrange(concat)
func.tempdiff <- as.matrix(dist(func.explanatory$avg.temp))
rownames(func.tempdiff) <- rownames(func.polmatrix)
colnames(func.tempdiff) <- rownames(func.polmatrix)
func.doydiff <- as.matrix(dist(func.explanatory$doy))
rownames(func.doydiff) <- rownames(func.polmatrix)
colnames(func.doydiff) <- rownames(func.polmatrix)
func.rounddiff <- as.matrix(dist(func.explanatory$round))
rownames(func.rounddiff) <- rownames(func.polmatrix)
colnames(func.rounddiff) <- rownames(func.polmatrix)
func.plantmatrix <- func.explanatory %>%
dplyr::select(site.date, concat) %>%
left_join(plants.sum, by = "site.date") %>%
dplyr::select(!site.date) %>%
pivot_wider(names_from = plant,
values_from = total_flowers,
values_fill = 0) %>%
arrange(concat) %>%
column_to_rownames(var = "concat")
func.plantbeta <- as.matrix(vegdist(func.plantmatrix, method = "horn"))
## Permutation Controls
func.concatdf <- data.frame(row = 1:nrow(func.polmatrix),
concat = rownames(func.polmatrix)) %>%
separate_wider_delim(cols = concat, names = c("date", "site", "round"), delim = ".") %>%
mutate(year = strftime(date, format = "%Y"),
site.year = paste(site, year, sep = "."))
func.control <- how(Within(type = "free"),
plots = Plots(as.factor(func.concatdf$site.year)),
blocks = as.factor(func.concatdf$round))
## MRM
set.seed(15723)
func.mrm <- suppressWarnings(MRM.base(as.dist(func.polbeta) ~
as.dist(func.tempdiff) +
as.dist(func.plantbeta) +
as.dist(func.doydiff),
method = "logistic",
nperm = 1000,
permute.control = func.control))
func.mrm
## $coef
## b pval
## (Intercept) -0.324750229 0.001
## as.dist.func.tempdiff. 0.017162610 0.096
## as.dist.func.plantbeta. 0.557472037 0.031
## as.dist.func.doydiff. 0.008174757 0.001
##
## $dev
## resid.dev resid.df dev.pval
## 7494.812 20096.000 1.000
family.polmatrix <- final.family %>%
group_by(concat, insect_family) %>%
summarise(num.active = sum(int)) %>%
pivot_wider(names_from = insect_family,
values_from = num.active,
values_fill = 0) %>%
arrange(concat) %>%
column_to_rownames(var = "concat")
family.polbeta <- as.matrix(vegdist(family.polmatrix, method = "horn"))
family.explanatory <- data.frame(concat = rownames(family.polmatrix)) %>%
separate(concat, into = c("date", "site", "round"), sep = "[.]", remove = F) %>%
mutate(date = as.Date(date),
round = as.numeric(round),
site.date = paste(site, date, sep = ".")) %>%
left_join(samptemps, by = c("date", "site", "round")) %>%
mutate(doy = strftime(date, format = "%j"),
year = strftime(date, format = "%Y")) %>%
arrange(concat)
family.tempdiff <- as.matrix(dist(family.explanatory$avg.temp))
rownames(family.tempdiff) <- rownames(family.polmatrix)
colnames(family.tempdiff) <- rownames(family.polmatrix)
family.doydiff <- as.matrix(dist(family.explanatory$doy))
rownames(family.doydiff) <- rownames(family.polmatrix)
colnames(family.doydiff) <- rownames(family.polmatrix)
family.rounddiff <- as.matrix(dist(family.explanatory$round))
rownames(family.rounddiff) <- rownames(family.polmatrix)
colnames(family.rounddiff) <- rownames(family.polmatrix)
family.plantmatrix <- family.explanatory %>%
dplyr::select(site.date, concat) %>%
left_join(plants.sum, by = "site.date") %>%
dplyr::select(!site.date) %>%
pivot_wider(names_from = plant,
values_from = total_flowers,
values_fill = 0) %>%
arrange(concat) %>%
column_to_rownames(var = "concat")
family.plantbeta <- as.matrix(vegdist(family.plantmatrix, method = "horn"))
## Permutation Controls
family.concatdf <- data.frame(row = 1:nrow(family.polmatrix),
concat = rownames(family.polmatrix)) %>%
separate_wider_delim(cols = concat, names = c("date", "site", "round"), delim = ".") %>%
mutate(year = strftime(date, format = "%Y"),
site.year = paste(site, year, sep = "."))
family.control <- how(Within(type = "free"),
plots = Plots(as.factor(family.concatdf$site.year)),
blocks = as.factor(family.concatdf$round))
## MRM
set.seed(71523)
family.mrm <- suppressWarnings(MRM.base(as.dist(family.polbeta) ~
as.dist(family.tempdiff) +
as.dist(family.plantbeta) +
as.dist(family.doydiff),
method = "logistic",
nperm = 1000,
permute.control = family.control))
family.mrm
## $coef
## b pval
## (Intercept) 0.097492788 1.000
## as.dist.family.tempdiff. 0.013756718 0.275
## as.dist.family.plantbeta. 0.506663633 0.034
## as.dist.family.doydiff. 0.008344169 0.001
##
## $dev
## resid.dev resid.df dev.pval
## 9335.089 19896.000 1.000
genus.polmatrix <- final.genus %>%
group_by(concat, insect_genus) %>%
summarise(num.active = sum(int)) %>%
pivot_wider(names_from = insect_genus,
values_from = num.active,
values_fill = 0) %>%
arrange(concat) %>%
column_to_rownames(var = "concat")
genus.polbeta <- as.matrix(vegdist(genus.polmatrix, method = "horn"))
genus.explanatory <- data.frame(concat = rownames(genus.polmatrix)) %>%
separate(concat, into = c("date", "site", "round"), sep = "[.]", remove = F) %>%
mutate(date = as.Date(date),
round = as.numeric(round),
site.date = paste(site, date, sep = ".")) %>%
left_join(samptemps, by = c("date", "site", "round")) %>%
mutate(doy = strftime(date, format = "%j"),
year = strftime(date, format = "%Y")) %>%
arrange(concat)
genus.tempdiff <- as.matrix(dist(genus.explanatory$avg.temp))
rownames(genus.tempdiff) <- rownames(genus.polmatrix)
colnames(genus.tempdiff) <- rownames(genus.polmatrix)
genus.doydiff <- as.matrix(dist(genus.explanatory$doy))
rownames(genus.doydiff) <- rownames(genus.polmatrix)
colnames(genus.doydiff) <- rownames(genus.polmatrix)
genus.rounddiff <- as.matrix(dist(genus.explanatory$round))
rownames(genus.rounddiff) <- rownames(genus.polmatrix)
colnames(genus.rounddiff) <- rownames(genus.polmatrix)
genus.plantmatrix <- genus.explanatory %>%
dplyr::select(site.date, concat) %>%
left_join(plants.sum, by = "site.date") %>%
dplyr::select(!site.date) %>%
pivot_wider(names_from = plant,
values_from = total_flowers,
values_fill = 0) %>%
arrange(concat) %>%
column_to_rownames(var = "concat")
genus.plantbeta <- as.matrix(vegdist(genus.plantmatrix, method = "horn"))
## Permutation Controls
genus.concatdf <- data.frame(row = 1:nrow(genus.polmatrix),
concat = rownames(genus.polmatrix)) %>%
separate_wider_delim(cols = concat, names = c("date", "site", "round"), delim = ".") %>%
mutate(year = strftime(date, format = "%Y"),
site.year = paste(site, year, sep = "."))
genus.control <- how(Within(type = "free"),
plots = Plots(as.factor(genus.concatdf$site.year)),
blocks = as.factor(genus.concatdf$round))
## MRM
set.seed(62223)
genus.mrm <- suppressWarnings(MRM.base(as.dist(genus.polbeta) ~
as.dist(genus.tempdiff) +
as.dist(genus.plantbeta) +
as.dist(genus.doydiff),
method = "logistic",
nperm = 1000,
permute.control = genus.control))
genus.mrm
## $coef
## b pval
## (Intercept) 0.203625449 1.000
## as.dist.genus.tempdiff. -0.008774277 0.518
## as.dist.genus.plantbeta. 0.627047758 0.018
## as.dist.genus.doydiff. 0.015189146 0.001
##
## $dev
## resid.dev resid.df dev.pval
## 11623.12 18911.00 1.00
We will check whether temperature affects the attractiveness of plants during our sampling periods. We will run this at all three levels of pollinator taxonomy again, just in case there are any differences associated with the different subsets of data.
# Functional Group
func.plants <- final.func %>%
mutate(round = as.factor(round)) %>%
group_by(date, site, round, avg.temp, plant) %>%
summarise(num.active = sum(int)) %>%
filter(num.active > 0)
totalplants <- func.plants %>%
group_by(date, site, round) %>%
summarise(total.active = sum(num.active))
plants.subset <- func.plants %>%
group_by(date, site, plant) %>%
tally() %>%
filter(n == 3) %>%
inner_join(func.plants, by = c("date", "site", "plant")) %>%
dplyr::select(-n) %>%
left_join(totalplants, by = c("date", "site", "round"))
activityplants <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|plant),
data = plants.subset,
family = binomial)
summary(activityplants)
## Family: binomial ( logit )
## Formula:
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,
## 2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |
## plant)
## Data: plants.subset
##
## AIC BIC logLik deviance df.resid
## 1544.2 1581.8 -762.1 1524.2 308
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## date (Intercept) 1.022156 1.01102
## site (Intercept) 0.267738 0.51743
## plant (Intercept) 6.595667 2.56820
## avg.temp 0.005845 0.07645 -0.94
## Number of obs: 318, groups: date, 47; site, 6; plant, 30
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.779903 1.888670 0.942 0.346
## poly(avg.temp, 2, raw = TRUE)1 -0.176164 0.130978 -1.345 0.179
## poly(avg.temp, 2, raw = TRUE)2 0.003366 0.002356 1.429 0.153
## round2 -0.004230 0.102757 -0.041 0.967
## round3 -0.043848 0.116240 -0.377 0.706
plants.resid <- simulateResiduals(activityplants, n = 250)
plot(plants.resid)
table1 <- final.genus %>%
group_by(site, year, date, round) %>%
tally() %>%
group_by(site, year) %>%
tally(name = "rounds.sampled") %>%
right_join(final.genus, by = c("site", "year"), multiple = "all") %>%
group_by(site, year, rounds.sampled) %>%
summarise(days.sampled = length(unique(date)))
table1
## # A tibble: 10 × 4
## # Groups: site, year [10]
## site year rounds.sampled days.sampled
## <chr> <chr> <int> <int>
## 1 Brush Creek Hill 2021 16 6
## 2 Brush Creek Hill 2022 15 6
## 3 Kebler Pass 2021 20 7
## 4 Kebler Pass 2022 18 7
## 5 Ohio Pass High 2022 19 7
## 6 Ohio Pass Low 2022 20 8
## 7 Rustlers Gulch 2021 21 7
## 8 Rustlers Gulch 2022 23 8
## 9 Stupid Falls 2021 20 7
## 10 Stupid Falls 2022 23 8
Relationship between round and temperature.
fig1 <- ggplot(samptemps %>% mutate(site.year = paste(site, year, sep = ".")),
aes(x = site.year, y = avg.temp, fill = as.factor(round))) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(color = "grey20",
alpha = 0.25,
show.legend = F,
position = position_jitterdodge(dodge.width = 0.75,
jitter.width = 0.2)) +
scale_fill_manual(values = brewer.pal(n = 8, "Dark2")[c(3, 5, 6)]) +
scale_color_manual(values = brewer.pal(n = 8, "Dark2")[c(3, 5, 6)]) +
labs(y = "Average Temperature (C)",
fill = "Round") +
theme(axis.text.x=element_text(angle = 60, hjust = 1),
text = element_text(size = 20),
axis.title.x = element_blank())
fig1
# ggsave("round-temp.png", fig1, width = 10, height = 8, units = "in")
General relationship between temperature and interaction occurrence.
First let’s calculate some broad-scale trends in interaction occurrence. Specifically, we will find the top 10 most frequently observed pollinator genera and the top 10 most frequently visited floral partners.
top.pols <- final.genus %>%
group_by(insect_genus) %>%
summarise(total.int = sum(int)) %>%
slice_max(order_by = total.int, n = 10)
top.plants <- final.genus %>%
group_by(plant) %>%
summarise(total.int = sum(int)) %>%
slice_max(order_by = total.int, n = 10)
top.fams <- final.family %>%
group_by(insect_family) %>%
summarise(total.int = sum(int)) %>%
slice_max(order_by = total.int, n = 10)
top.func <- final.func %>%
group_by(gross_ID) %>%
summarise(total.int = sum(int)) %>%
slice_max(order_by = total.int, n = 6)
top.ints <- final.genus %>%
filter(insect_genus %in% top.pols$insect_genus,
plant %in% top.plants$plant) %>%
group_by(insect_genus, plant) %>%
summarise(total.int = sum(int)) %>%
filter(total.int >= 5)
top.ints2 <- final.func %>%
group_by(int.ID) %>%
summarise(total.int = sum(int)) %>%
slice_max(order_by = total.int, n = 50)
# Bombus
bombus.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
filter(group == "2")
bombus.df[bombus.df$plant == "Potentilla pulcherrima",]$x[which.max(bombus.df[bombus.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 28.506
bombus.df[bombus.df$plant == "Heliomeris multiflora",]$x[which.max(bombus.df[bombus.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 35.918
bombus.df[bombus.df$plant == "Senecio crassulus",]$x[which.max(bombus.df[bombus.df$plant == "Senecio crassulus",]$predicted)]
## [1] 32.188
bombus.df[bombus.df$plant == "Erigeron speciosus",]$x[which.max(bombus.df[bombus.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 33.221
bombus.plot <- ggplot(bombus.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_manual(values = brewer.pal(n = 5, "Dark2")[c(1, 2, 3, 5)]) +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
coord_cartesian(ylim = c(0, 0.004)) +
theme(text = element_text(size = 20),
legend.position = "none")
bombus.plot
ggsave("bombus-predcurves.png", bombus.plot, width = 6, height = 6, units = "in")
# Syrphidae
syrphid.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
filter(group == "2")
syrphid.df[syrphid.df$plant == "Potentilla pulcherrima",]$x[which.max(syrphid.df[syrphid.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 26.543
syrphid.df[syrphid.df$plant == "Heliomeris multiflora",]$x[which.max(syrphid.df[syrphid.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 33.914
syrphid.df[syrphid.df$plant == "Senecio crassulus",]$x[which.max(syrphid.df[syrphid.df$plant == "Senecio crassulus",]$predicted)]
## [1] 30.014
syrphid.df[syrphid.df$plant == "Erigeron speciosus",]$x[which.max(syrphid.df[syrphid.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 30.965
syrphid.df[syrphid.df$plant == "Pseudocymopterus montanus",]$x[which.max(syrphid.df[syrphid.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 24.545
syrphid.plot <- ggplot(syrphid.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_brewer(palette = "Dark2") +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
coord_cartesian(ylim = c(0, 0.004)) +
theme(text = element_text(size = 20),
legend.position = "none")
syrphid.plot
ggsave("syrphid-predcurves.png", syrphid.plot, width = 6, height = 6, units = "in")
# SBSB
sbsb.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
filter(group == "2")
sbsb.df[sbsb.df$plant == "Potentilla pulcherrima",]$x[which.max(sbsb.df[sbsb.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 18.806
sbsb.df[sbsb.df$plant == "Heliomeris multiflora",]$x[which.max(sbsb.df[sbsb.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 26
sbsb.df[sbsb.df$plant == "Senecio crassulus",]$x[which.max(sbsb.df[sbsb.df$plant == "Senecio crassulus",]$predicted)]
## [1] 22.238
sbsb.df[sbsb.df$plant == "Erigeron speciosus",]$x[which.max(sbsb.df[sbsb.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 23.198
sbsb.df[sbsb.df$plant == "Pseudocymopterus montanus",]$x[which.max(sbsb.df[sbsb.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 17.089
sbsb.plot <- ggplot(sbsb.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_brewer(palette = "Dark2") +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
coord_cartesian(ylim = c(0, 0.004)) +
theme(text = element_text(size = 20),
legend.position = "none")
sbsb.plot
ggsave("sbsb-predcurves.png", sbsb.plot, width = 6, height = 6, units = "in")
# Other Fly
fly.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
filter(group == "2")
fly.df[fly.df$plant == "Potentilla pulcherrima",]$x[which.max(fly.df[fly.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 19.758
fly.df[fly.df$plant == "Heliomeris multiflora",]$x[which.max(fly.df[fly.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 27.084
fly.df[fly.df$plant == "Senecio crassulus",]$x[which.max(fly.df[fly.df$plant == "Senecio crassulus",]$predicted)]
## [1] 23.34
fly.df[fly.df$plant == "Erigeron speciosus",]$x[which.max(fly.df[fly.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 24.306
fly.df[fly.df$plant == "Pseudocymopterus montanus",]$x[which.max(fly.df[fly.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 17.855
fly.plot <- ggplot(fly.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_brewer(palette = "Dark2") +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
coord_cartesian(ylim = c(0, 0.004)) +
theme(text = element_text(size = 20),
legend.position = "none")
fly.plot
ggsave("fly-predcurves.png", fly.plot, width = 6, height = 6, units = "in")
# Beetles
beetles.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
filter(group == "2")
beetles.df[beetles.df$plant == "Potentilla pulcherrima",]$x[which.max(beetles.df[beetles.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 22.525
beetles.df[beetles.df$plant == "Senecio crassulus",]$x[which.max(beetles.df[beetles.df$plant == "Senecio crassulus",]$predicted)]
## [1] 26
beetles.df[beetles.df$plant == "Erigeron speciosus",]$x[which.max(beetles.df[beetles.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 27.084
beetles.df[beetles.df$plant == "Pseudocymopterus montanus",]$x[which.max(beetles.df[beetles.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 20.424
beetles.plot <- ggplot(beetles.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_manual(values = brewer.pal(n = 5, "Dark2")[c(1, 3, 4, 5)]) +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
coord_cartesian(ylim = c(0, 0.004)) +
theme(text = element_text(size = 20),
legend.position = "none")
beetles.plot
ggsave("beetles-predcurves.png", beetles.plot, width = 6, height = 6, units = "in")
# Hemipteran
hemip.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
filter(group == "2")
hemip.df[hemip.df$plant == "Potentilla pulcherrima",]$x[which.max(hemip.df[hemip.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 25.806
hemip.df[hemip.df$plant == "Heliomeris multiflora",]$x[which.max(hemip.df[hemip.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 33.225
hemip.df[hemip.df$plant == "Senecio crassulus",]$x[which.max(hemip.df[hemip.df$plant == "Senecio crassulus",]$predicted)]
## [1] 29.356
hemip.df[hemip.df$plant == "Erigeron speciosus",]$x[which.max(hemip.df[hemip.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 30.361
hemip.df[hemip.df$plant == "Pseudocymopterus montanus",]$x[which.max(hemip.df[hemip.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 23.919
hemip.plot <- ggplot(hemip.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_brewer(palette = "Dark2") +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
coord_cartesian(ylim = c(0, 0.004)) +
theme(text = element_text(size = 20),
legend.position = "none")
hemip.plot
ggsave("hemip-predcurves.png", hemip.plot, width = 6, height = 6, units = "in")
legend.plot <- ggplot(hemip.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
geom_line(linewidth = 1) +
scale_color_brewer(palette = "Dark2") +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability",
color = "Floral Partner",
linetype = "Floral Partner") +
coord_cartesian(ylim = c(0, 0.011)) +
theme(text = element_text(size = 18),
legend.position = "bottom") +
guides(color = guide_legend(nrow = 2),
linetype = guide_legend(nrow = 2))
ggsave("predcurves-legend.png", legend.plot, width = 10, height = 6, units = "in")
# SBSB
## Halictidae
halictidae.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.family, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(insect_family = "Halictidae", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
filter(group == "2")
fam.halic <- final.family %>%
filter(insect_family == "Halictidae") %>%
uncount(weights = int)
halictidae.df[halictidae.df$plant == "Potentilla pulcherrima",]$x[which.max(halictidae.df[halictidae.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 17.089
mcplot.halictidae <- ggplot(halictidae.df,
aes(x = x)) +
geom_line(aes(y = predicted),
color = "#7570B3",
linewidth = 1.5) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
fill = "grey",
color = "grey",
alpha = 0.1,
linetype = "dotted") +
geom_vline(xintercept = 17.089,
color = "black",
linetype = "dashed") +
geom_rug(data = fam.halic,
aes(x = avg.temp),
sides = "b") +
coord_cartesian(xlim = c(15, 38),
ylim = c(0, 0.006)) +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability",
color = "Plant",
fill = "Plant") +
theme(text = element_text(size = 20),
legend.position = "bottom")
mcplot.halictidae
# ggsave("mcplot-halictidae.png", mcplot.halictidae, width = 6, height = 6, units = "in")
## Andrenidae
andrenidae.df <- data.frame() %>%
bind_rows(bind_cols(ggpredict(mc.family, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(insect_family = "Andrenidae", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
filter(group == "2")
fam.andren <- final.family %>%
filter(insect_family == "Andrenidae") %>%
uncount(weights = int)
andrenidae.df[andrenidae.df$plant == "Potentilla pulcherrima",]$x[which.max(andrenidae.df[andrenidae.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 19.948
mcplot.andrenidae <- ggplot(andrenidae.df, aes(x = x, )) +
geom_line(aes(y = predicted),
color = "#7570B3",
linewidth = 1.5) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
fill = "grey",
color = "grey",
alpha = 0.1,
linetype = "dotted") +
geom_vline(xintercept = 19.948,
color = "black",
linetype = "dashed") +
geom_rug(data = fam.andren,
aes(x = avg.temp),
sides = "b") +
coord_cartesian(xlim = c(15, 38),
ylim = c(0, 0.006)) +
labs(x = "Temperature (C)",
y = "Predicted Interaction Probability") +
theme(text = element_text(size = 20))
mcplot.andrenidae
# ggsave("mcplot-andrenidae.png", mcplot.andrenidae, width = 6, height = 6, units = "in")
# Kebler Pass, 2021-07-21
kp.072121.bombus <- final.genus %>%
filter(insect_genus == "Bombus",
site == "Kebler Pass",
date == "2021-07-21") %>%
group_by(site, date, round, insect_genus) %>%
summarise(genus.int = sum(int)) %>%
left_join(final.genus, by = c("site", "date", "round", "insect_genus")) %>%
group_by(round, insect_genus, plant, avg.temp, genus.int) %>%
summarise(int = sum(int),
rel.int = int/genus.int) %>%
filter(rel.int > 0)
kp.072121.plant <- plants %>%
group_by(date, site) %>%
summarise(total_site = sum(total_flowers)) %>%
left_join(plants, by = c("date", "site")) %>%
filter(site == "Kebler Pass",
date == "2021-07-21") %>%
group_by(site, plant) %>%
summarise(rel.abund = total_flowers/total_site) %>%
mutate(plant = replace(plant, plant %in% c("Agoseris glauca", "Boechera stricta", "Collomia linearis", "Linum lewisii", "Helianthella quinquenervis", "Potentilla pulcherrima"), "0ther"))
kp.072121.plot <- ggplot(kp.072121.bombus,
aes(x = avg.temp, y = rel.int, fill = plant)) +
geom_bar(position = "stack",
stat = "identity",
col = "grey50",
width = 0.6) +
scale_x_continuous(breaks=c(24.0, 30.7, 31.6)) +
scale_fill_manual(values = brewer.pal(n = 9, "Blues")[c(3, 6, 8, 9)]) +
labs(x = "Sampling Temperature (C)",
y = "Interaction Occurrence",
col = "Plant") +
theme(text = element_text(size = 20),
legend.position = "bottom",
axis.title = element_blank(),
legend.title = element_blank()) +
coord_cartesian(xlim = c(22, 32))
kp.072121.plot
kp.072121.plantplot <- ggplot(kp.072121.plant,
aes(x = site, y = rel.abund, fill = plant)) +
geom_bar(position = "stack",
stat = "identity",
width = 1.2) +
scale_fill_manual(labels = c("Other", sort(unique(kp.072121.plant$plant))[2:5]),
values = c("lightgrey", "#C6DBEF", "#4292C6", "#08519C", "#08306B"))
# Rustlers Gulch, 2021-07-27
rg.072721.bombus <- final.genus %>%
filter(insect_genus == "Bombus",
site == "Rustlers Gulch",
date == "2021-07-27") %>%
group_by(site, date, round, insect_genus) %>%
summarise(genus.int = sum(int)) %>%
left_join(final.genus, by = c("site", "date", "round", "insect_genus")) %>%
group_by(round, insect_genus, plant, avg.temp, genus.int) %>%
summarise(int = sum(int),
rel.int = int/genus.int) %>%
filter(rel.int > 0)
rg.072721.plant <- plants %>%
group_by(date, site) %>%
summarise(total_site = sum(total_flowers)) %>%
left_join(plants, by = c("date", "site")) %>%
filter(site == "Rustlers Gulch",
date == "2021-07-27") %>%
group_by(site, plant) %>%
summarise(rel.abund = total_flowers/total_site) %>%
mutate(plant = replace(plant, plant %in% c("Artemisia dracunculus", "Erigeron elatior", "Potentilla fruticosa", "Pseudocymopterus montanus", "Senecio bigelovii", "Valeriana edulis"), "0ther"))
rg.072721.plot <- ggplot(rg.072721.bombus %>%
filter(insect_genus == "Bombus"),
aes(x = avg.temp, y = rel.int, fill = plant)) +
geom_bar(position = "stack",
stat = "identity",
col = "grey50",
width = 0.6) +
scale_x_continuous(breaks=c(22.5, 25.1, 29.1)) +
scale_fill_manual(values = brewer.pal(n = 9, "Blues")[c(1, 2, 5, 7)]) +
labs(x = "Sampling Temperature (C)",
y = "Interaction Occurrence",
col = "Plant") +
theme(text = element_text(size = 20),
legend.position = "bottom",
axis.title = element_blank(),
legend.title = element_blank()) +
coord_cartesian(xlim = c(22, 32))
rg.072721.plot
rg.072721.plantplot <- ggplot(rg.072721.plant,
aes(x = site, y = rel.abund, fill = plant)) +
geom_bar(position = "stack",
stat = "identity",
width = 1.2) +
scale_fill_manual(labels = c("Other", sort(unique(rg.072721.plant$plant))[2:6]),
values = c("lightgrey", "#F7FBFF", "#DEEBF7", "#C6DBEF", "#6BAED6", "#2171B5"))